﻿namespace GeostatisticalTool.Classes
{
    using System;
    using System.Collections;
    using DotSpatial.Data;
    using DotSpatial.Topology;
    using DotSpatial.Projections;
    using DotSpatial.Controls;
    using System.Collections.Generic;


    class LocalTrend : AbstractInterpolator, IInterpolator
    {

        private int order = 0;
        private string[] serie = null;
        private Matriz result = null;
        private double meanValue = 0;

        public LocalTrend(IFeatureSet shapeLayer, int idField, int order) :
            base(shapeLayer, idField)
        {
            this.method = "LocalTrend";
            this.order = order;
            if (order >= 1 || order <= 3)
            {
                serie = SerieCoef(order);
            }
        }

        public LocalTrend(IFeatureSet shapeLayer, int idField, int order, double meanValue) :
            base(shapeLayer, idField)
        {
            this.method = "LocalTrend";
            this.meanValue = meanValue;
            this.order = order;
            if (order >= 1 || order <= 3)
            {
                serie = SerieCoef(order);
              //  GenerateMatrix(order, this.npoints, this.Kpoints);
            }
        }


        new public ParameterInterpolator GetParameters()
        {
            ParameterInterpolator par = new ParameterInterpolator();
            par.Index = this.order;
            return par;
        }

        new public void SetParameters(ParameterInterpolator parameters)
        {
            this.order = parameters.Index;
            if (order >= 1 || order <= 3)
            {
                serie = SerieCoef(order);
                GenerateMatrix(order, this.npoints, this.Kpoints);
            }

        }



        public double Seq(int numPoints,
            int order, int col, int row, double x, double y)
        {

            string coe1 = serie[col];
            string coe2 = serie[row];

            //if (coe1 == "0" && coe2 == "0")
            //{ return numPoints; }

            double val1 = Oper(coe1, x, y);
            double val2 = Oper(coe2, x, y);

            return val1 * val2;


        }

        public string SeqM(int numPoints,
            int order, int col, int row, double x, double y)
        {

            string coe1 = serie[col];
            string coe2 = serie[row];

            //if (coe1 == "0" && coe2 == "0")
            //{ return "c"; }


            return coe1 + "*" + coe2;


        }


        public double Seq(int numPoints,
          int order, int col, int row, double x, double y, double z)
        {

            string coe1 = serie[col];

            //if (coe1 == "0")
            //{ return z; }

            double val1 = Oper(coe1, x, y);

            return val1 * z;


        }


        public double Oper(
            string type, double val1, double val2)
        {
            switch (type)
            {
                case "x":
                    return val1;
                case "y":
                    return val2;
                case "xy":
                    return val1 * val2;
                case "x2":
                    return val1 * val1;
                case "y2":
                    return val2 * val2;
                case "x2y":
                    return val1 * val1 * val2;
                case "xy2":
                    return val2 * val2 * val1;
                case "x3":
                    return val1 * val1 * val1;
                case "y3":
                    return val2 * val2 * val2;
                default:
                    return 1.0;
            }
        }

        public string[] SerieCoef(
            int order)
        {
            if (order == 1)
                return new string[3] { "0", "x", "y" };
            if (order == 2)
                return new string[6] { "0", "x", "y", "x2", "xy", "y2" };
            if (order == 3)
                return new string[10] { "0", "x", "y", "x2", "xy", "y2", "x3", "x2y", "xy2", "y3" };

            return new string[3] { "0", "x", "y" };
        }


        public void GenerateMatrix(int order, int num, double[,] data)
        {
            int matrixSize = SerieCoef(order).Length;
            Matriz matrix = new Matriz(matrixSize);
            Matriz ZMatrix = new Matriz(matrixSize, 1);
            double width = this.extent[2] - this.extent[0];
            double height = this.extent[3] - this.extent[1];
            double difZ = this.extentZ[1] - this.extentZ[0];
            for (int v = 0; v < num; v++)
            {
                double x = (data[v, 0] - extent[0]) / width;
                double y = (data[v, 1] - extent[1]) / height;
                double z = data[v, 2];
                for (int row = 0; row < matrixSize; row++)
                {
                    for (int col = 0; col < matrixSize; col++)
                    {
                        matrix[row + 1, col + 1] = matrix[row + 1, col + 1] + Seq(num, order, col, row, x, y);
                    }
                    ZMatrix[row + 1, 1] = ZMatrix[row + 1, 1] + Seq(num, order, row, 0, x, y, z);
                }
            }

            Matriz inver = matrix.Inverse();
            Matriz ind = matrix * inver;

            result = inver * ZMatrix;
        }


        public void GenerateMatrix(int order, int num, List<Kpoint> data)
        {
            int matrixSize = SerieCoef(order).Length;
            Matriz matrix = new Matriz(matrixSize);
            Matriz ZMatrix = new Matriz(matrixSize, 1);
            double width = this.extent[2] - this.extent[0];
            double height = this.extent[3] - this.extent[1];
            double difZ = this.extentZ[1] - this.extentZ[0];
            for (int v = 0; v < num; v++)
            {
                double x = (data[v].X - extent[0]) / width;
                double y = (data[v].Y - extent[1]) / height;
                double z = data[v].W;
                for (int row = 0; row < matrixSize; row++)
                {
                    for (int col = 0; col < matrixSize; col++)
                    {
                        matrix[row + 1, col + 1] = matrix[row + 1, col + 1] + Seq(num, order, col, row, x, y);
                    }
                    ZMatrix[row + 1, 1] = ZMatrix[row + 1, 1] + Seq(num, order, row, 0, x, y, z);
                }
            }

            Matriz inver = matrix.Inverse();
            Matriz ind = matrix * inver;

            result = inver * ZMatrix;
        }


        public Matriz GenerateMatrixCross(int order, int num, double[,] data)
        {
            int matrixSize = SerieCoef(order).Length;
            Matriz matrix = new Matriz(matrixSize);
            Matriz ZMatrix = new Matriz(matrixSize, 1);
            double width = extent[2] - extent[0];
            double height = extent[3] - extent[1];
            double difZ = extentZ[1] - extentZ[0];
            for (int v = 0; v < num; v++)
            {
                double x = (data[v, 0] - extent[0]) / width;
                double y = (data[v, 1] - extent[1]) / height;
                double z = data[v, 2];
                for (int row = 0; row < matrixSize; row++)
                {
                    for (int col = 0; col < matrixSize; col++)
                    {
                        matrix[row + 1, col + 1] = matrix[row + 1, col + 1] + Seq(num, order, col, row, x, y);
                    }
                    ZMatrix[row + 1, 1] = ZMatrix[row + 1, 1] + Seq(num, order, row, 0, x, y, z);
                }
            }

            Matriz inver = matrix.Inverse();
            Matriz ind = matrix * inver;

            Matriz result = inver * ZMatrix;
            return result;
        }

        public Matriz GenerateMatrixCross(int order, int num, List<Kpoint> data)
        {
            int matrixSize = SerieCoef(order).Length;
            Matriz matrix = new Matriz(matrixSize);
            Matriz ZMatrix = new Matriz(matrixSize, 1);
            double width = extent[2] - extent[0];
            double height = extent[3] - extent[1];
            double difZ = extentZ[1] - extentZ[0];
            for (int v = 0; v < num; v++)
            {
                double x = (data[v].X - extent[0]) / width;
                double y = (data[v].Y - extent[1]) / height;
                double z = data[v].W;
                for (int row = 0; row < matrixSize; row++)
                {
                    for (int col = 0; col < matrixSize; col++)
                    {
                        matrix[row + 1, col + 1] = matrix[row + 1, col + 1] + Seq(num, order, col, row, x, y);
                    }
                    ZMatrix[row + 1, 1] = ZMatrix[row + 1, 1] + Seq(num, order, row, 0, x, y, z);
                }
            }

            Matriz inver = matrix.Inverse();
            Matriz ind = matrix * inver;

            Matriz result = inver * ZMatrix;
            return result;
        }

        public override double[] Interpolate(double x, double y, double z, bool include)
        {
            if (this.order == 0)
            {
                return new double[1] { this.meanValue };
            }

            List<Kpoint> cpointsF;

            this.Search.ListPointAvailables(out cpointsF, x, y,z, include);

            if (cpointsF.Count == 1)
            {
                return new double[2] { cpointsF[0].W, double.NaN };
            }

            if (cpointsF.Count == 0)
            {
                return new double[2] { double.NaN, double.NaN };
            }



            GenerateMatrix(order, cpointsF.Count, cpointsF);

            double sum = 0;
            int i = 0;
            double width = extent[2] - extent[0];
            double height = extent[3] - extent[1];
            double difZ = extentZ[1] - extentZ[0];

            foreach (string s in serie)
            {
                double valOpe = Oper(s, (x - extent[0]) / width, (y - extent[1]) / height);
                double resPara = 0;
                resPara = result[i + 1, 1].Realnumber;
                sum = sum + (resPara * valOpe);
                i++;
            }
            double val = (sum * difZ) + extentZ[0];
            return new double[1] { sum };
        }

        public override double[] CrossValidation()
        {

            double[] output = new double[2];
            List<Kpoint> cpointsCross = new List<Kpoint>();
            double SSerror = 0;
            double sumAllValues = 0;
            double quad = 0; ;

            for (int drop = 0; drop < this.NumPoints; drop++)
            {
                cpointsCross = new List<Kpoint>();
                for (int i = 0; i < this.NumPoints; i++)
                    if (i != drop)
                        cpointsCross.Add(this.Kpoints[i]);

               // this.result = GenerateMatrixCross(this.order, this.npoints - 1, cpointsCross);
                Kpoint value = this.Kpoints[drop];
                double[] inter = this.Interpolate(value.X, value.Y,value.Z, false);
                value.Zinte = inter[0];
                value.SdtError = value.W - value.Zinte;
                sumAllValues += value.SdtError;
                quad += (value.Zinte * value.Zinte);
                SSerror += (value.SdtError * value.SdtError);

            }

            output[0] = Math.Sqrt(SSerror / this.npoints); //RMS
            output[1] = sumAllValues / this.npoints; //Mean Prediction errors
            this.GenerateMatrix(this.order, this.npoints, this.Kpoints);

            return output;
        }


    }
}